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13a  abstract  (Maximum  2GC  wo/xs) 

We  use  direct  numerical  simulations  to  analyze  the  evolution  of  a  temporally 
growing  two-dimensional  free  shear  layer  seeded  with  a  dilute  suspension  of 
bubbles  under  gravity.  The  bubble  concentrations  sure  dilute  enough  so  tlmt  bubble- 
bubble  interactions  Eire  negligible,  but  cumulative  effects  of  bubbles  alter  the  flow 
field.  The  evolution  of  the  bubbles  can  then  be  determined  by  tracking  many 
individual  bubbles,  and  the  flow  field  is  advanced  using  the  Navier-Stokes 
equations  with  a  coupling  term  in  the  momentum  equation  representing  the  effect 
of  the  bubbles  on  the  flow.  We  interpret  the  results  in  terms  of  the  difference  in  the 
vorticity,  bubble  concentration,  and  pressure  fields  relative  to  the  passive  or  one¬ 
way  coupled  case.  Due  to  the  nature  of  the  vorticity  production  mechanism,  the 
net  circulation  is  not  affected  by  the  bubbles,  but  local  variations  do  occur, 
especially  near  the  vortex  center.  In  addition  to  the  effect  of  the  bubbles  on  the 
flow,  the  bubble  field  is  also  altered  as  a  result  of  the  two-way  coupling.  The 
‘  location  of  bubble  accumulation  is  shifted  away  from  the  vortex  center  and  the 
magnitude  of  this  accumulation  is  reduced  relative  to  the  passive  case. 


1A.  SUBJECT  TERMS 


bubbly  flow,  dilute  concentration,  two-way  coupling, 
Navier-Stokes  equations,  numerioal  simulation 


17.  SECURITY  CLASSIfKATION 
Of  REPORT 

Unclassified 


NSN  7540-01-2a0-5500 


IS.  SECURITY  CLASSIPICATION 
OP  THIS  PACE 

Unclassified 


19.  SECURITY  CLASSIPICATION 
OP  ABSTRACT 

Unclassified 


15.  NUMBER  CP  PAGES 


16.  PRICE  CODE 


20.  UMITATION  OP  A3STRAC 


Stancara  roem  298  (3*v.  2-59) 

Dw  APEtJ  tia. 


Two-way  Coupling  in  Shear  Layers  with 
Dilute  Bubble  Concentrations 


G.  R.  Ruetsch*  and  E.  Meiburg 


Department  of  Aerospace  Engineering, 

University  of  Southern  California,  Los  Angeles,  CA  90089-1191 


“Present  aiddress:  Center  for  Turbulence  Research, 
Stanford  University,  Stanford,  CA  94305 


Sixbaaitted^o  Physics  of  Fluids  A,  September  1993 


Abstract 


Accesioh  For 

N'TIS  CRA&I 
OTIC  TAD 

U. ..i.-ifv.ni 

Ci':.::!: 


By . 


Ju.I  I  / 


Dist 


/I- 1 


Avail  ci' 


Spec 


We  use  direct  numerical  simulations  to  analyze  the  evolution  of  a  temporally  grow¬ 
ing  two-dimensional  shear  layer  seeded  with  dilute  concentrations  of  bubbles  under 
gravity.  The  bubble  concentrations  are  dilute  enough  so  that  bubble-bubble  inter¬ 
actions  can  be  neglected,  but  are  large  enough  fur  cumulative  effects  of  bubbles  to 
influence  the  flow.  The  evolution  of  the  bubble  field  is  determined  by  tracking  many 
individual  bubbles,  and  the  flow  field  is  advanced  by  using  the  Navier-Stokes  equations 
with  a  coupling  term  representing  the  effect  of  the  bubbles  on  the  flow.  We  inter¬ 
pret  the  results  in  terms  of  the  vorticity,  density,  and  pressure  fields  relative  to  the 
one-way  coupled  or  passive  case.  For  the  coupled  case  we  observe  a  reduction  in  the 
magnitude  of  the  vorticity  and  pressure  gradients  near  the  vortex  center.  In  addition 
to  modification  of  the  flow,  we  observe  that  the  accumulation  of  bubbles  is  smaller 
and  the  location  of  the  equilibrium  points  are  shifted  further  from  the  vortex  center 
as  a  result  of  the  coupling.  We  explore  how  these  changes  are  modified  by  different 
Froude  numbers  and  bubble  sizes.  The  differences  between  passive  and  coupled  cases 
usually  increase  due  to  larger  accumrdations  as  we  consider  larger  bubbles.  However, 
for  certain  Froude  numbers  an  optimum  coupling  is  observed  at  intermediate  bubble 
sizes  due  to  the  absence  of  equilibrium  points  for  large  bubbles. 


PACS  numbers:  47.55.Kf,  47.25.Jn 
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1  Introduction 


One  important  and  common  feature  which  occurs  in  a  variety  of  fluid  systems  is  the  pres¬ 
ence  of  shear  layers.  The  evolution  of  shear  layers  has  received  much  attention  by  researchers 
due  to  the  critical  role  it  plays  in  mixing,  momentum  transport,  amd  transition  to  turbu¬ 
lence.  As  a  result,  there  have  been  many  advances  in  understamding  the  nature  of  single 
phase  shear  layers,  such  as  the  observation  of  coherent  spanwise  structures  which  arise  from 
Kelvin-Helmholtz  instabilities,^’^  and  the  three-dimensional  structures  which  arise  when 
these  spanwise  structures  become  unstable.^’^ 

Recently,  the  behavior  of  small  particles  in  shear  layers  has  been  inv^tigated.  The 
applications  of  such  studies  are  far  reaching,  ranging  &om  the  mixing  of  fuel  in  engines  to 
the  dispersion  of  pollutants  in  the  environment.  The  nature  of  this  advectiou  depends  on 
the  density  of  the  particles  relative  to  the  fluid  of  the  shear  layer.  On  one  hand,  dispersion 
of  heavy  or  aerosol  particles  occurs  when  these  particles  are  introduced  into  shear  layers, 
which  is  observed  both  in  experiments®’®  and  simulations.^’^  On  the  other  hand,  particles 
which  are  lighter  than  the  carrier  fluid,  such  as  air  bubbles  in  water,  accumiflate  near  the 
center  of  vortices  in  cellular  flows^  and  in  temporally  evolving  shear  layers.^®  One  common 
assumption  made  in  all  the  previously  mentioned  numerical  studies  is  to  only  consider  dilute 
suspensions  of  pairticles.  This  assumption  allows  one  to  neglect  the  interactions  between 
particles  and  the  fluid  and  among  particles  themselves. 

There  have  been  several  approaches  to  studying  non-dilute  particle  laden  flows.  For 
example,  Biesheuvel  and  Gorrisen^^  use  a  kinetic  theory  approach  to  derive  one-dimensional 
conservation  equations  in  their  investigation  of  void  fraction  distiurbances  for  the  case  of 
large  bubbles  at  large  void  fractions.  Cook  and  Harlow^^  employed  ensemble-averaged  two- 
field  equations  with  various  closure  models  to  investigate  a  bubble-laden  von  Karman  vortex 
street. 

In  this  study,  we  wish  to  consider  “weakly-dilute”  suspensions  of  bubbles,  namely  flows 
with  bubbles  that  are  dilute  enough  so  that  bubble-bubble  interactions  can  be  neglected, 
but  not  dilute  enough  to  ignore  cumulative  effects  of  bubbles  on  the  flow.  We  also  wish  to 
do  this  in  a  direct  fashion,  to  avoid  making  assumptions  about  closure  models.  A  similar  ap¬ 
proach  has  been  used  for  heavy  particles  in  the  simiilations  of  Squires  and  Eaton^^  and  most 
recently  Elghobashi  and  Truesdell,^^  who  investigated  turbulence  modification  by  particles 
by  including  a  source  of  momentum^ in  the  Navier- Stokes  equations  which  accounts  for  the 
net  force  on  the  particles  back  on  the  fluid.  We  use  this  approach  for  the  case  of  bubbles, 
which  due  to  their  massless  nature  yield  a  different  contribution  to  the  momentiim  equations 
of  the  flow. 

In  the  next  section,  we  discuss  the  governing  equations  for  the  bubble  motion  and  then 
derive  the  momentum  equations  for  the  flow  making  use  of  the  weakly-dilute  aissumption 
in  a  generad  context.  We  then  discuss  the  particulzir  flow  configuration  we  pursue  in  this 
study.  A  brief  discussion  of  the  numerical  methods  and  methods  for  analysis  of  the  data  are 
presented,  followed  by  the  results  of  the  simulations. 
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2  Governing  equations 

The  flow  we  are  considering  concerns  the  motion  of  a  dilute  concentration  of  gas  bubbles  in 
a  liquid,  such  as  water.  Our  approach  to  this  problem  is  to  solve  equations  governing  each 
phase  which  incorporate  the  effects  from  the  other  phase,  resulting  in  a  two-way  coupling. 
Throughout  this  section  we  refer  to  the  gas  bubbles  as  either  the  bubble  or  particle  phase 
and  the  water  or  liquid  as  the  fluid  phase. 

Unlike  other  studies,  our  goal  is  not  to  derive  field  equations  for  both  phases.  We  derive 
field  equations  only  for  the  fluid  phase,  and  calculate  the  bubble  phase  by  following  discrete 
bubbles.  The  equations  discussed  in  this  section  apply  to  any  weakly-dilute  bubbly  flow.  We 
specify  these  equations  to  our  particular  flow  field  and  provide  initial  conditions  in  a  later 
section. 


2.1  Bubble  Equation 

We  begin  by  discussing  the  equations  used  to  calculate  the  bubble  trajectories  in  the  fluid. 
Depending  on  the  flow  conditions  around  the  bubble,  there  are  several  different  equations  one 
can  use  here.  In  this  study  we  consider  bubbles  small  enough  such  that  they  are  dominated 
by  viscous  forces,  and  consider  surface  tension  large  enough  to  assume  the  bubbles  spherical. 
In  addition,  we  assiime  that  the  bubbles  maintain  a  constant  volume,  so  that  we  can  consider 
these  bubbles  as  rigid  spheres.  Maxey  and  Riley have  derived  such  an  equation  for  particles 
in  general  which,  neglecting  the  Faxen  corrections  for  a  non-uniform  flow,  is  given  by: 

Duj  dVj 

Dt  dt 

-|-67ro/i  (u,-  —  Vi)  +  67ro^/i  /  (1) 

-'o  -  r) 


dVi  (  \  X  .  ^  / 


where  Ui{x j,t)  and  Vi(t)  are  the  fluid  and  pau'ticle  velocities,  m/  amd  are  the  mass  of 
the  fluid  displaced  by  the  particle  and  the  mass  of  the  particle,  o  is  the  particle  radius, 
and  fi  is  the  fluid  viscosity.  Here  we  have  assumed  that  the  acceleration  of  gravity  g  points 
in  the  negative  13— direction.  The  terms  on  the  right-hand  side  of  EJq.  (1)  represent  the 
gravitational  force,  pressure  force  in  absence  of  the  particle,  the  added  mass  effects  of  the 
form  given  by  Auton  et  Stokes  drag,  and  the  Basset  history  term.  In  addition  to 
neglecting  the  Faxen  corrections,  we  also  choose  to  neglect  the  Basset  history  term.  The 
basis  for  doing  so  lies  in  the  assumptions  used  to  derive  Eq.  (1),  where  the  Reynolds  number 
based  on  the  particle  radius  and  slip  velocity  is  zero.  Departures  from  this  condition  have 
been  determined  to  result  in  a  more  quickly  decaying  kernel  of  the  Basset  history  term.^”^ 
We  nondimensionalize  Eq.  (1)  with  U  and  5  (not  to  be  confused  with  the  Kronecker  delta, 
6ij,  used  above  in  the  gravitational  term)  as  the  velocity  and  length  scales,  and  using  the 
notation  of  Maxey^  we  obtain  the  following  equation  for  the  acceleration  of  the  particle; 

f  =  ^(u..  +  W«.-V5)  +  |k^  (2) 
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We  have  introduced  three  nondimensional  parameters  in  Eq.  (2),  A,  and  W.  The  settling 
velocity  parameter,  W,  is  defined  as; 

W  =  (”*/  ~  ^p)9 
6i:afiU 

and  represents  the  ratio  of  gravitational  to  viscous  effects,  where  WU  gives  the  terminal 
rise/settling  velocity  of  a  particle  in  a  stiU  fluid.  The  mass  ratio  parameter,  Ti,  is  defined  as: 


h  + 

and  reflects  the  differences  in  the  fluid  and  particle  densities.  The  material  derivative  in 
Eq.  (2)  represents  the  effect  of  the  presstire  gradients  on  the  particle  motion,  and  therefore 
■R  plays  a  crucial  role  in  determining  how  particle  and  fluid  element  trajectories  differ.  The 
mass  ratio  peirzuneter  can  cover  the  range  of  0  <  7^  <  2.  For  “R,  =  2/3,  the  pressure  has  the 
same  effect  on  the  particle  as  on  a  fluid  element,  and  corresponds  to  the  case  of  neutrally 
buoyant  particles.  For  7^  =  0,  corresponding  to  heavy  particles,  the  pressure  forces  have  no 
effect  on  the  particle  motion.  For  bubbles  with  7^  =  2,  the  pressure  forces  are  three  times  as 
important  to  the  bubble  motion  relative  to  a  fluid  element.  It  is  for  this  reason  that  heavy 
particles  are  dispersed  by  vortices  while  bubbles  are  entrapped  by  vortices.  In  this  study,  we 
consider  only  bubbles  with  72.  =  2. 

The  inertia  parameter,  A,  is  defined  as; 

^ _ 6irafi6 

(m,  +  |m/)  U 


and  is  the  inverse  of  the  Stokes  number.  The  inertia  parameter  reflects  the  relative  impor¬ 
tance  of  viscous  to  inertial  effects,  md  in  this  study  is  large  reflecting  the  dominance  of  the 
viscous  forces. 

At  this  point  we  should  comment  on  the  credibility  of  using  the  Stokes  drag  law  in 
Eqs.  (1)  and  (2).  One  requisite  for  using  this  expression  is  that  the  Reynolds  number  based 
on  the  slip  velocity  and  particle  dieuneter  be  small: 

V 

In  absence  of  gravitational  effects  this  restriction  is  met  for  the  values  of  A  we  consider  here, 
but  since  we  are  interested  in  cases  where  W  ^  0  this  condition  can  be  violated.  We  can 
correct  this  using  an  empirical  coefficient  based  on  Re,  from  Clift  tt 


f  =  1.0  + 0.15Rel^^ 


which  results  in  the  following  equation: 

^  =  +  +  (3) 

In  addition  to  altering  the  drag  term,  the  non-zero  rise  velocity  further  justifies  the  elimina¬ 
tion  of  the  Basset-history  term  from  this  equation  due  to  finite  Reynolds  number  effects. 

In  developing  Eq.  (3)  we  have  assumed  that  only  a  single  bubble  is  present.  But  since  we 
are  dealing  with  a  weakly-dilute  suspension  of  bubbles,  the  bubble-bubble  interactions  are 
neglected,  and  we  can  calculate  each  bubble  trajectory  using  Eq.  (3)  independently. 
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2.2  Fluid  Equations 


In  this  section  we  derive  the  governing  equations  which  describe  the  evolution  of  the  fluid 
phase.  Since  we  are  concerned  with  the  motion  of  air  bubbles  in  an  incompressible  liquid, 
we  begin  the  anzdysis  of  the  governing  equations  for  the  fluid  phase  with  the  incompressible 
Navier-Stokes  equations: 

dui 


Pi- 


Dui 

W 


dP 


d^Ui 


=  -5^-W<3i  +  (‘-5;r  +  ‘f 


(5) 


where  P  is  the  pressure,  n  is  the  dynamic  viscosity,  and  is  the  force  of  the  bubbles  on  the 
fluid.  We  can  determine  bi  from  from  the  equation  of  motion  for  a  single  bubble,  Eq.  (1), 
which  can  be  restated  as: 

=  Vi  +  5,-  =  0 


where  the  last  equality  results  from  substituting  mp  =  0  in  the  first  term.  Here  Vi  and  Si 
represent  the  net  body  and  surface  forces  on  the  bubble.  From  Eq.  (1)  we  obtain  Vi  =  mjgSzi, 
and  Si  consists  of  adl  the  other  terms  on  the  right-hand  side  of  the  equation.  It  is  important 
to  differentiate  between  these  body  amd  surface  forces,  since  Newton’s  third  law  of  action- 
reaction  implies  that  the  force  exerted  by  a  single  bubble  on  the  fluid  is  simply 


Bi  =  -Si  =  Vi 

=  PfCgSzi 

where  C  =  |iro^  is  the  volume  of  the  bubble.  The  coupling  term  bi  in  the  momentum 
equation  exists  at  the  bubble  interface,  but  rather  than  calculate  the  interaction  between 
phcises  in  such  a  fashion  we  choose  a  more  macroscopic  approach.  For  N  bubbles  imder  the 
we£Lkly-dilute  assumption  we  obtain  the  overaJl  force  as 

NpfCgSzi. 

where  C  is  the  average  volume  of  the  bubbles.  This  extensive  qumtity  can  then  be  converted 
to  the  average  intensive  force  bi  by  dividing  by  the  total  volume: 

-  hi  =  fjCpfgS^i 

where  rj  is  the  global  or  average  number  density,  defined  as  the  number  of  bubbles  per  volume. 
We  C2in  also  use  the  average  void  fraction, 

€  =  ijC 


in  this  expression: 


bi 


^Pjghi 


Although  this  force  applies  only  at  the  fluid-bubble  interface,  for  omi  macroscopic  approach 
we  can  define  the  local  force  of  the  bubbles  on  the  fluid  as: 


hi  —  tpjgSzi 
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where  6  is  the  local  void  haction.  Using  this  expression  for  the  force  per  unit  volume  of  the 
bubbles  on  the  fluid  we  can  write  £q.  (5)  as: 

Dui  dP  d^Ui 

'’'oT ="5;: 

We  can  interpret  the  last  term  of  Eq.  (6)  in  terms  of  the  density  of  the  fluid-bubble  mixtxire, 
which  is  given  by: 

p  =  pf{l  -  e)  (7) 

From  Eqs.  (6)  and  (7)  we  see  that  the  net  effect  of  the  weakly-dilute  suspension  of  bubbles 
is  to  introduce  a  buoyancy  term  which  is  a  function  of  the  local  void  fraction. 

The  next  step  is  to  determine  the  reference  state  of  the  fluid  obtained  by  setting  u,-  =  0 
amd  e  =  e  in  Eq.  (6): 

For  cases  with  c  =  0  and  thus  bi  =  0,  we  obtain  the  familiar  hydrostatic  equation,  where 
Po  —  pj.  However,  in  our  case  this  procedure  3rields  a  reference  density  given  by: 

Po  =  P/{1  -  e) 

or  simply  the  average  density  of  the  fluid-bubble  mixture.  Using  this  relation  with  Eq.  (7) 
we  can  define  the  perturbation  density  of  the  mixture  as: 

/  =  P-Po 

=  /»/(?-  e)  (8) 

If  we  now  make  the  substitution  P  =  Po+P',  subtract  the  reference  state  from  Eq.  (6),  divide 
by  />/,  and  non-dimensionalize  Eqs.  (4),  (6),  and  (8)  according  to  the  following  variables: 

u^  =  ulU',  x'^  =  xl5]  —tUfS 
p^  =  pIpj\  =  €/c;  =  PliPfU^) 

we  obtain  the  following  set  of  equations,  with  the  ‘•"’s  omitted; 


|2i  =  o 

azi 

(9) 

Dui 

dP'  ,  1  d^ui  1  ,, 

(10) 

W~' 

dxi'^Redxj 

P'  =  e(l  -  c) 

(11) 

In  addition  to  including  e  in  this  set  of  equations,  we  have  also  introduced  two  additional 
parameters,  the  Reynolds  and  Froude  numbers,  given  by: 

Rc  =  U8lv-  T=UI{gS)'^l'^ 
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We  should  point  out  several  features  concerning  Eqs.  (9)-(ll).  The  most  obvious  is  that 
Eq.  (10)  is  the  same  equation  that  one  obtains  in  a  fluid  of  variable  density  under  the 
Boussinesq  approximation.  The  Boussinesq  approximation  assumes  that  density  variations 
are  small  enough  so  that  the  density  appears  as  a  constant  in  all  terms  except  the  buoyancy 
term.  These  are  the  same  assumptions  we  made  when  specifying  our  fundamental  equations, 
so  the  fact  that  we  obtain  the  Boussinesq  equations  is  no  surprise.  The  difference  between 
a  Boussinesq  fluid  juad  our  two-phase  flow  is  that  in  the  case  of  the  former  the  density  or 
temperature  field  is  determined  by  a  convection-diffusion  equation,  whereas  for  the  two-phase 
flow  the  density  field  is  a  function  of  the  void  fraction,  Eq.  (11),  which  is  in  turn  dependent 
on  the  equation  of  motion  for  the  bubbles,  Eq.  (3). 

3  Flow  Configuration 

In  the  previous  section  we  discussed  the  set  of  equations  needed  to  solve  for  each  component 
of  a  two-phase  flow.  In  developing  these  equations  we  left  the  flow  length  md  velocity  scales, 
8  and  C/,  undetermined,  since  we  wished  to  emphasize  the  fact  that  these  equations  are  valid 
for  any  weakly-dilute  bubbly  flow.  Having  done  this,  we  now  turn  our  attention  to  the 
specific  flow  we  consider. 

The  flow  configuration  we  study  is  a  two-dimensional  shear  layer  initially  with  a  parallel 
flow  in  the  ii— direction  and  gravity  acting  in  the  negative  I3— direction.  This  base  flow 
has  a  tanh  velocity  profile,  and  is  perturbed  using  eigenfunctions  corresponding  to  the  most 
unstable  mode  of  inviscid  theory.^^  The  velocity  scale  U  is  the  velocity  difference  across  the 
shear  layer,  and  the  length  scale  8  is  taken  to  be  the  vorticity  thickness  defined  as: 

_ I _ . 

For  these  scales,  the  nondimensionad  initial  profile  is  given  by: 

u(x3)  =  ^  tanh(2z3) 

and  the  initiad  perturbation  is  given  with  the  fundamental  wave  number  a  =  0.8892.  We 
give  this  perturbation  an  amplitude  of  0.01. 

In  addition  to  defining  the  flow'scales,  we  must  also  discuss  our  selection  of  the  non- 
dimensional  parameters  appearing  in  the  finad  set  of  equations  we  use  to  solve  the  two-phase 
flow.  This  final  set  of  equations  consists  of  Eq.  (3)  for  the  bubble  phase  and  Eqs.  (9)-(ll) 
for  the  fluid-bubble  mixture.  From  these  equations  we  have  i?e,  and  c  appearing  in  the 
fluid  equations  and  A,  W,  and  Tl  appearing  in  the  bubble  equation  as  our  non-dimensional 
parameters.  The  mass  ratio  parameter  is  fixed  at  “R  :=  2,  and  to  simplify  matters  we 
consider  an  overall  void  fraction  of  c  =  0.01  md  a  Reynolds  number  of  Rt  =  1000  in  all 
cases.  This  leaves  us  with  A,  and  W.  However,  these  three  parameters  cannot  be 
chosen  independently.  We  should  recognize  that  both  A  and  Re  give  ratios  of  viscous  to 
inertial  forces,  likewise  W  (actually  WA)  and  ^  relate  gravitational  to  inertial  effects.  The 
difference  between  these  two  sets  of  parameters  concerns  the  length  scales  on  which  these 
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terms  operate.  Therefore,  if  we  include  the  ratio  of  bubble  to  flow  length  scales,  ajS,  in  our 
list  of  parameters,  we  can  then  express  the  bubble  parauneters  A  and  W  in  terms  of  the  flow 
parameters  Re.  and  T‘. 

A^ 

-=!(?)’ a-')  I 

As  a  resTilt  of  these  relations,  for  fixed  6,  Rt,  and  R  we  have  only  two  independent  parameters, 
and  we  choose  these  to  be  T  and  either  ajh  or  A.  By  choosing  a/5  or  -4  as  a  global  paranneter 
we  au:e  restricting  ourselves  to  simulations  with  uniform  bubble  sizes.  However,  we  could 
easily  assign  individual  values  to  a/5  (also  requiring  different  values  for  A  and  W)  for  each 
bubble  in  our  simulations. 


(12) 


4  Numericad  simulation 


The  direct  numerical  simulation  used  to  calculate  the  evolution  of  the  flow  is  not  unlike 
methods  previously  used  in  Rayleigh-Benard  convection,^®  or  more  precisely  the  strati¬ 
fied  sheair-layer  code  of  Wang  and  Maxey.^^  We  use  a  pseudospectral  technique  to  advamce 
Eqs.  (9) — (11)  in  our  two-dimensional  simulations  on  a  64  x  128  grid.  We  simulate  a  tem¬ 
poral  shear  layer,  and  thus  apply  periodic  boundary  conditions  in  the  ii-direction.  Because 
we  consider  bubbles  with  a  finite  rise  velocity  that  can  affect  the  flow,  the  no  normal  flow 
conditions  at  the  vertical  boundjunes  commonly  used  in  the  Rayleigh-Benard  simulations  are 
not  applicable  here.  Instead,  we  apply  periodic  boundary  conditions  in  the  13— direction. 
We  do  this  by  extending  the  domain  in  X3  so  that  flow  in  the  center  of  the  simulation  is  not 
influenced  by  the  bovmdeiry,  and  then  apply  an  additional  shear  layer  of  opposite  sense  to  the 
main  shear  layer  at  the  boundaries  to  allow  for  periodicity.  The  shear  layer  at  the  boundary 
has  twice  the  thickness  of  the  real  shear  layer,  which  in  addition  to  the  decay  of  the  initial 
perturbation  away  from  the  central  shear  layer  stabilizes  the  layer  across  the  boundary.  The 
box  size  is  then  £  =  27r/Q  ~  7.1  in  the  zj— direction  and  2L  in  the  Z3— direction. 

We  rewrite  Eq.  (10)  as  follows: 


dui 

dt 


=  (ux«),-^ 


1  d^Ui 


(14) 


and  advance  the  flow  field  using  a  second-order  Adams-Beishforth  scheme  on  the  non-linear 
and  buoymcy  terms  amd  a  second-order  Crank-Nicolson  scheme  on  the  linear  viscous  terms. 
The  pressure  term  and  the  continuity  condition  of  Eq.  (9)  are  taken  caze  of  by  projecting 
the  Fourier  coefficients  of  the  velocity  field  onto  an  incompressible  space: 


Ui  = 


where  ki  is  the  wavenumber  vector. 

We  initially  seed  the  entire  flow  field  uniformly  with  bubbles  of  the  same  size,  so  that  the 
spacing  between  bubbles  in  both  Zi—  and  Z3— directions  is  the  same.  The  bubble  locations 
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are  axivanced  applying  the  second  order  predictor-corrector  scheme  used  in  Ruetsch  and 
Meiburg^^  to  Eq.  (3)  after  the  flow  field  is  advanced.  This  scheme  requires  only  a  fixst- 
order  accurate  estimate  of  the  material  derivative  of  the  fluid  velocity  at  both  predictor  and 
corrector  steps  in  order  to  obtain  an  overall  second-order  scheme.  In  order  to  determine  the 
fluid  velocity  at  the  bubble  locations,  the  Hermite  interpolation  scheme  of  Balachandar  and 
Maxey  is  used. 

The  remmning  computational  issue  concerns  the  method  for  obtaining  the  values  for  the 
density  at  the  grid  points  from  the  bubble  locations.  For  simulations  with  bubbles  of  uniform 
volume,  such  as  the  simulations  we  consider  in  this  study,  we  have  the  relation: 


e  = 


where  7;  has  been  nondimensionalized  by  and  thus  the  perturbation  density  in  Eq.  (11) 
can  be  expressed  as 

o'  =  c(l  -  77), 

so  that  we  only  need  to  establish  a  method  of  determining  the  number  density  at  each  grid 
point.  In  order  to  do  so  efficiently,  we  use  an  interpolation  scheme  where  a  bubble  only 
affects  the  four  grid  points  that  define  the  cell  in  which  the  bubble  is  located.  In  general,  the 
contribution  of  each  bubble  to  these  grid  points  is  weighted  according  to  the  location  of  the 
bubble  within  the  grid  ceil  and  the  volume  of  the  bubble,  but  for  our  uniform  bubble  size, 
this  contribution  is  weighted  only  by  the  location.  For  simplicity  we  use  a  linear  weighting 
in  each  direction.  Thus,  if  a  bubble  is  located  in  a  grid  cell  with  the  bottom  left  comer 
referenced  by  the  grid  coordinates  (i,  j),  with  the  fractional  distance  of  the  bubble  from 
this  grid  point  (A*,,  A*,),  then  the  contributions  from  this  bubble  to  the  non-dimensional 
number  density  is  given  as  follows: 


77, j  (1  —  An)(l  -  ^xz)NaRlDlNjOT 

Vi+i.j  ^*1(1  -  ^x3)Ngrid I Ntot 

GRID  jN JOT 

ViJ+l  (1  - 


where  Ngrid  =  64  x  128  is  the  total  number  of  grid  points  and  Njoj  is  the  total  number  of 
bubbles  in  the  flow.  For  an  equally  spaced  grid,  we  cam  think  of  the  number  density  as  the 
number  of  bubbles/grid  cell,  in  which  caise  the  average  number  density  is  given  by  the  ratio 
Njot INGRID,  and  the  local  number  density  is  determined  by  summing  the  (1  — An)(l  — Arj), 
etc.,  terms.  The  ratio  of  these  gives  the  nondimensional  number  density.  Even  for  large 
numbers  of  bubbles  this  techniques  creates  high  frequencies  in  the  Fourier  coefficients  of  the 
density  field,  and  therefore  requires  some  sort  of  filtering.  We  apply  an  exponential  filter  to 
the  Fourier  coefficients  of  the  density  field: 


where  7  is  chosen  so  exp( — 7)  gives  the  machine  accuracy,  and  kmax  is  the  maximum 
wavenumber  in  either  direction. 
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Although  we  simulate  a  two-dimensional  flow,  the  void  fraction  here  represents  a  volu¬ 
metric  void  fraction  which  can  be  interpreted  by  assuming  the  bubble  field  is  periodic  in  the 
xj— direction  with  the  same  initial  spacing  between  bubbles  in  X2  as  in  ii  and  X3,  with  the 
flow  uniform  in  the  Xj— direction.  We  would  like  to  choose  a  large  number  of  bubbles  which 
result  in  smoother  density  profiles,  but  at  the  same  time  require  that  £  be  small.  Because  we 
assume  that  there  is  no  interaction  between  bubbles  in  our  weakly-dilute  cases,  we  satisfy 
both  conditions  by  carrying  “virtual”  bubbles  along  with  the  “actual”  bubbles.  As  a  result 
we  need  to  differentiate  between  the  actual  number  of  bubbles,  N act ^  and  the  total  number 
Ntot  —  ^ACT  +  Nvifix-  We  can  then  obtain  a  smoother  density  values  at  the  grid  points 
by  summing  over  all  the  bubbles  and  normaJizing  by  Ntot,  however  we  determine  the  other 
parameters  of  the  flow  based  on  Nact-  As  an  example,  we  determine  the  bubble  radius  from 
£  and  Nact-  In  turn,  the  bubble  parameters  A  and  W  are  given  by  this  bubble  radius  and 
Eqs.  (12)  and  (13).  A  detailed  summary  of  the  parameters  for  all  simulations  is  given  in 
Table  I. 

5  Analysis 


We  begin  our  analysis  of  the  effect  of  bubbles  on  the  shear  layer  by  investigating  the  changes 
which  occur  to  the  vorticity  field.  The  general  vorticity  equation  for  a  Boussinesq  fluid  with 
gravity  acting  in  the  negative  13-direction  is  given  by: 


dui  1 


1  d^Ui 
Re  dxj 


The  first  two  terms  on  the  right-hand  side  of  this  equation  are  vorticity  production  terms 
resulting  from  vortex  stretching/tilting  and  horizontal  density  fluctuations  being  acted  on 
by  gravity.  The  last  term  represents  the  viscous  diffusion  of  vorticity.  In  our  case,  where  we 
are  dealing  with  only  two  spatial  dimensions,  the  vorticity  equation  reduces  to: 

Dt  ~  Redx]  ^  ^ 

where  u  here  is  in  the  ij— direction.  Therefore,  for  a  two-dimensional  flow,  the  only  way 
the  vorticity  of  a  material  element  can  change  is  either  by  diffusion  or  due  to  coupling  term 
as  a  result  of  the  horizontal  density  gradients.  Because  the  density  field  is  periodic  in  the 
xi-direction,  the  net  production  of  vorticity  across  the  central  shear  layer  resulting  from  the 
coupling  term  is  zero.  Therefore,  diffusion  is  the  only  process  which  effects  the  circulation 
across  the  shear  layer. 

In  addition  to  the  effect  of  the  bubbles  on  the  shear  layer,  the  two-way  coupling  will  also 
affect  aspects  of  the  bubble  motion  relative  to  the  passive  case.  Since  the  deviation  of  bubble 
motion  from  fluid  particles  is  dominated  by  the  effects  of  pressure,  such  as  the  accumulation 
of  bubbles  into  regions  of  low  pressure,  the  effects  of  coupling  on  the  bubble  motion  can 
be  exeimined  from  the  equation  governing  the  pressure  field.  This  equation  is  obtained  by 
taking  the  divergence  of  Eq.  (14): 


d  ,  ,  1  a/ 
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where,  in  addition  to  the  nonlinear  terms,  the  density  perturbation  provides  an  extra  source 
to  the  Poisson  equation  for  the  pressure.  It  is  interesting  to  note  from  Eqs.  (15)  and  (16) 
that  the  horizontal  gradients  in  the  density  affect  the  vorticity,  while  it  is  the  vertical  density 
gradients  that  Jiffect  the  pressure. 

6  Results 

Throughout  this  section  we  axe  interested  not  so  much  in  describing  how  the  flow  and  bubble 
motion  evolve  with  two-way  coupling  in  absolute  terms,  but  rather  how  this  evolution  differs 
from  the  one-way  coupled  or  passive  case.  The  idea  here  is  to  provide  some  insight  into  the 
differences  between  these  two  cases  which  can  be  used  to  develop  models  that  predict  the 
coupled  results  based  on  the  knowledge  of  the  passive  state.  It  is  for  this  re<ison  that  when  we 
discuss  results  from  coupled  simulations,  we  always  compare  these  to  simulations  performed 
with  passive  bubbles.  The  passive  simulations  are  nm  by  setting  ^  =  oo  which  eliminates 
the  effect  of  the  bubbles  on  the  flow,  but  by  keeping  all  other  pauameters,  including  W, 
identical  to  the  coupled  case.  We  then  can  generate  spatial  fields  and  statistics  based  on  the 
differences  between  the  two  cases. 

Before  we  discuss  the  results  from  the  coupled  simulations  and  how  they  differ  from  their 
passive  counterparts,  we  must  first  review  some  basic  featmes  of  what  occurs  in  the  passive 
case.  Although  the  density  fields  for  the  passive  cases  will  vary  depending  on  A  and  W,  the 
evolution  of  the  vorticity  field  will  not  change.  We  show  a  time  sequence  of  the  vorticity 
field  for  the  passive  case  in  Fig.  1.  From  this  sequence  we  observe  that  the  flow  remains 
roughly  parallel  up  to  t  ~  10,  at  which  time  we  observe  the  formation  of  the  vortex  core  due 
to  the  Kelvin- Helmholtz  instability.  During  the  latter  stages  of  the  shear  layer  development 
we  observe  a  roughly  circular  nature  of  the  flow  near  the  vortex  center.  Simple  model  flows 
reflecting  this  circular  nature  have  been  used  to  investigate  the  motion  of  bubbles  in  vortical 
flows.  One  such  model  which  allows  a  simple  analytical  solution  to  bubble  trajectories  is 
that  of  a  solid-body  vortex.^^  Two  important  featxires  obtained  from  this  solution  are  the 
equilibrium  points,  or  locations  where  the  bubbles  come  to  rest,  and  the  rate  at  which  the 
bubbles  aire  captured  by  the  equilibrium  points.  Without  any  drag  correction,  an  exponential 
rate  of  entrapment  or  accumulation  is  observed.  This  entrapment  is  a  function  of  A  alone, 
eind  reaches  an  optimal  vcilue  at  ~  1.  In  the  present  study  we  consider  ^  1,  and 
therefore  the  general  trend  we  expect  is  for  the  accumulation  of  bubbles  to  be  greater  as  A 
becomes  smaller.  In  the  previous  work^®  we  differentiated  between  the  accumulation  rate  of 
bubbles  near  the  vortex  center  and  the  ability  of  the  vortex  to  trap  bubbles  for  cases  with 
rise  velocities.  The  accumulation  rate  near  the  vortex  center,  as  we  mentioned  previously, 
obtains  a  maodmum  at  A  ~  1,  whereas  we  find  the  that  the  ability  of  the  vortex  to  trap 
bubbles  throughout  the  fluid  increases  with  decreasing  A.  However,  for  the  range  of  A  we 
consider  here,  when  we  discuss  entrapment  or  accumulation  we  axe  referring  to  both  these 
phenomena. 

The  location  of  the  equilibrium  points  is  a  fimction  of  both  A  and  W.  For  W  =  0  the 
equilibrium  point  coincides  with  the  vortex  center.  As  W  increases,  the  equilibrium  point 
moves  away  from  the  vortex  center  along  a  line  whose  slope  depends  on  A.  For  large  A,  this 
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line  is  roughly  horizontal,  corresponding  to  a  force  b2dance  at  the  equilibrium  point  primarily 
between  the  drag  and  gravitational  forces.  For  smaller  A,  the  line  becomes  more  vertical  as 
pressure  forces  become  more  important  in  the  force  balance. 

When  the  slip  velocity  of  the  bubble  becomes  large,  which  is  most  noticeable  for  cases  with 
rise  velocities,  we  must  apply  corrections  to  the  Stokes  drag  assumption  as  in  Eq.  (3).  In  such 
cases,  the  analytical  solution  for  the  bubble  trajectories  discussed  above  is  no  longer  valid. 
One  important  difference  between  the  anal)rtical  solution  for  cases  without  drag  correction 
and  the  actual  bubble  trajectories  is  the  change  in  entrapment  rate  ne2ir  the  equilibrium 
point.  The  entrapment  rate  becomes  a  fimction  of  both  A  and  VV,  with  the  entrapment  rate 
decreasing  as  VV  increases.  This  reduced  accumulation  rate  can  be  understood  by  noting 
that  the  coeflBicient  of  the  drag  term  in  Eq.  (3)  is  fA,  where  the  correction  /  is  a  function  of 
the  slip  velocity.  Because  the  equilibrium  point  moves  further  away  from  the  vortex  center 
with  increasing  W,  the  slip  velocity  at  the  equilibrium  point,  (=  u,),  and  hence  /  incre«ise 
accordingly.  Therefore,  when  correcting  for  the  deviation  from  Stokes  drag  we  are  effectively 
considering  an  inertia  parameter  of  Aactual  =  IbA,  where  fs  is  the  correction  at  the 
equilibrium  point.  Therefore,  we  expect  the  bubble  accumulation  to  vary  inversely  with 
both  A  and  W. 

Having  reviewed  these  characteristics  of  the  passive  bubble  motion  in  a  solid-body  vortex, 
we  now  proceed  to  discuss  the  results  of  the  two-way  coupled  simulations  of  a  bubbly  shear 
layer.  In  the  following  sections  we  group  or  results  according  to  different  Froude  numbers, 
varying  the  bubble  size  within  each  section. 

6.1  T  ^  1:  Intermediate  case 

We  begin  our  analysis  for  the  T  —  I  cases  by  considering  the  laurgest  bubble  size,  correspond¬ 
ing  to  an  inertia  parameter  of  =  5.8.  Up  to  about  t  =  20,  we  observe  small  differences 
between  the  passive  and  coupled  runs  since  during  this  time  the  vortex  core  and  pressure 
gradient  have  not  developed  to  a  point  where  the  bubble  accumulation  is  large  enough  to 
affect  the  flow.  At  t  =  30,  however,  we  do  observe  differences  between  the  two  cases.  A  plot 
of  the  vorticity,  density,  amd  pressure  for  both  passive  and  coupled  cases  is  shown  in  Fig.  2a. 
It  is  helpful  in  comparing  the  coupled  and  passive  cases  to  plot  the  difference  between  these 
two  fields  obtained  by  subtracting  the  two  fields  point  by  point.  These  differences  are  also 
shown  in  Fig.  2a.  Consistent  with  observations  from  the  paissive  bubble  motion  in  a  solid- 
body  vortex,  we  observe  the  accumulation  of  bubbles  (signified  by  a  large  negative  density) 
to  the  right  of  the  vortex  center.  We  adso  observe  that  this  accumulation  occurs  as  a  result  of 
the  depletion  of  bubbles  from  the  braid  region  of  the  flow,  or  region  between  the  vortex  cores. 
From  Eq.  (15)  we  know  that  such  accumulations  of  bubbles  result  in  vorticity  production 
on  either  side  of  the  acciimulation  given  by  f~^{dp'ldx{).  To  the  left  of  the  accumulation 
this  term  generates  negative  vorticity  and  to  the  right  positive  vorticity,  resulting  in  the 
observed  decrease  of  vorticity  in  the  vortex  center  and  the  increase  in  vorticity  along  the 
right-heind  edge  of  the  vortex.  Once  again  we  should  emphasize  that  the  net  effect  of  the 
production  term  is  zero,  so  that  the  global  circulation  is  preserved,  but  this  production  does 
create  localized  variations. 

At  later  times,  such  as  in  Fig.  2b,  we  observe  a  greater  change  in  the  vorticity  between  the 
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pjtssive  and  coupled  cases.  Similar  to  earlier  times,  the  vorticity  near  the  center  is  reduced 
while  the  vorticity  near  the  edge  of  the  vortex  is  increased.  Unlike  earlier  times,  we  note 
that  the  magnitude  of  this  difference  is  greatest  near  the  vortex  core,  and  a  weaker  but 
more  voluminous  (positive)  difference  surrounds  the  vortex.  Statistically,  this  difference  is 
observed  in  the  PDF  of  the  vorticity  difference,  Auj,  shown  in  Fig.  3.  In  this  figure,  where 
the  vorticity  difference  is  normalized  by  its  rms  value,  Au',  we  observe  in  addition  to  the 
large  contribution  to  the  PDF  from  small  vzdues  of  the  vorticity  difference  a  small  peak  of 

negative  vorticity  difference  occurring  at - 10 Aw'.  It  is  these  values  which  correspond 

to  the  large  vorticity  differences  near  the  center  of  the  vortex.  There  iire  two  reasons  for 
the  occurrence  of  intense  negative  vorticity  difference  and  the  absence  of  a  similar  region 
of  positive  vorticity  difference.  The  first  is  due  to  the  shape  of  the  accumiilation  region. 
This  region  hcis  a  sharper  gradient  on  the  left  side,  near  the  center  of  the  vortex,  than  on 
the  right  side.  Although  this  accumulation  region  produces  no  net  vorticity,  the  negative 
vorticity  resiilting  from  the  density  gradient  is  more  concentrated.  The  second  reason  for  the 
absence  of  a  large  positive  vorticity  difference  results  from  the  difference  in  the  locations  of 
the  equilibrium  point  of  the  bubble  field  and  the  stagnation  point  of  the  flow  field.  Although 
the  location  of  bubble  accumulation  remains  approximately  stationary  throughout  the  flow 
evolution,  the  production  of  vorticity  due  to  this  accumulation  applies  to  material  elements 
of  the  fluid.  The  stagnation  point  in  the  vortex  core  occurs  to  the  left  of  the  accumulation 
region,  and  therefore  the  negative  vorticity  production  occurring  here  is  applied  to  the  same 
fluid  elements  resulting  in  the  large  negative  vorticity  difference  neax  the  vortex  core.  To  the 
right  of  the  vortex  core  the  positive  vorticity  production  affects  different  fluid  elements  as 
they  are  swept  past  the  accrimulation.  This  sweeping  results  in  a  smaller  vorticity  difference 
applied  to  a  larger  area. 

In  addition  to  the  effect  of  the  bubbles  on  the  flow,  the  modification  of  the  flow  has  an 
effect  on  the  entrapment  of  bubbles.  This  is  more  clearly  seen  at  later  times,  for  example 
in  Fig.  2b.  The  difference  in  the  density  field  cam  be  attributed  to  the  reduction  in  the 
magnitude  of  the  pressure  at  the  vortex  center.  The  smaller  pressure  gradients  affect  both 
the  location  and  the  entrapment  rate  into  the  equilibrium  points.  The  entrapment  rate  is 
reduced,  amd  the  location  is  shifted  away  from  the  vortex  center. 

Up  to  this  point  we  have  discussed  the  basic  mechanisms  concerning  the  modification  of 
both  flow  and  bubble  fields  resulting  from  two-way  coupling  for  a  single  set  of  parameters. 
As  we  change  the  bubble  size,  and  in  later  sections  the  Froude  number,  we  observe  approxi¬ 
mately  the  same  behavior.  We  therefore  focus  more  on  the  quantitative  changes  rather  than 
qualitative  features  resulting  from  different  pairameters.  When  discussing  how  different  pa¬ 
rameters  tiffect  the  results,  we  must  first  discuss  what  effect  the  change  in  parameters  has  on 
the  passive  simulations,  amd  then  proceed  to  the  coupled  simulations  their  difference  relative 
to  the  passive  case.  Keeping  these  ideas  in  mind,  we  now  look  at  simulations  with  smaller 
bubbles. 

We  now  consider  smaller  bubbles  with  A  =  10.3  and  23.  In  order  to  keep  T  constant, 
this  increase  in  A  is  accompanied  by  a  decrezise  in  W.  For  the  case  of  passive  bubbles  in 
a  solid- body  vortex  we  know  that  the  effect  of  the  bubble  size  on  the  entrapment  rate  is 
two-fold.  On  one  hand  the  entrapment  rate  decreases  with  increasing  A  for  the  range  of  A 
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considered  here.  However,  we  also  observe  greater  accumulation  rates  when  W  decreases, 
since  the  effective  inertia  parameter  near  the  equilibrium  point  is  f^A  and  /e  decreases  with 
W.  In  addition  to  the  effect  of  bubble  size  on  the  entrapment  rate,  we  also  observe  a  shift  in 
the  location  of  the  equilibrium  point  towards  the  vortex  center  as  the  bubble  size  decreases. 
We  can  clearly  see  these  two  effects  in  Fig.  4  for  passive  bubbles.  The  accumulation  is  closer 
to  the  center  of  the  vortex  relative  than  the  passive  case  in  Fig.  2b.  However,  in  spite  of  the 
reduction  of  the  effective  inertia  parameter  due  to  the  smaller  drag  correction,  we  observe 
a  reduced  accumulation  rate,  apparent  &om  the  minimum  in  the  p*  field.  The  reduction  in 
drag  correction  does  not  offset  the  change  in  A. 

Similar  to  the  case  of  A  =  5.8,  we  observe  that  the  accumulation  is  smaller  and  occurs 
further  away  from  the  vortex  center  for  the  coupled  simulation  than  in  the  passive  simulation 
for  these  smaller  bubbles.  However,  the  degree  to  which  this  happens  decreases  with  smadler 
bubbles.  This  difference  is  quantified  in  Fig.  5,  which  shows  time  series  of  both  the  maximum 
and  minimum  number  densities  for  all  bubble  sizes  in  both  passive  and  coupled  simulations. 
(The  minimum  in  p'  corresponds  to  the  maximum  in  t;.)  From  this  figure  we  see  that  the 
accumulation  decreases  with  bubble  size  for  the  passive  and  coupled  cases.  In  addition, 
the  difference  between  the  passive  and  coupled  cases  also  decreases  with  bubble  size.  This 
is  what  one  expects,  the  largest  difference  between  passive  and  coupled  cases  occurs  when 
the  density  gradient  and  hence  accumulation  is  lairgest,  which  in  turn  occurs  for  the  larger 
bubbles.  Note  that  at  later  times  for  the  case  where  A  —  5.8  the  reduction  of  vorticity  is  so 
great  that  the  vortex  becomes  unable  to  hold  ail  the  bubbles  it  previously  entrapped,  and  as 
a  result  the  number  density  decreases,  as  shown  in  Fig.  5.  As  a  restilt,  the  density  minimum 
for  A  =  5.8  is  only  slightly  larger  than  for  A  —  10.3,  which  at  the  last  time  frames  show  a 
roughly  equivalent  maximum  vorticity  difference. 

Aside  from  the  2D  contour  plots  of  the  density  fields  and  the  time  series  of  the  maximum 
and  minimum  number  density,  another  useful  quantity  to  exeunine  is  the  vertical  profile  of  the 
density  field,  obtained  by  averaging  the  density  over  horizontal  slices.  The  advantage  of  these 
profiles  over  the  full  2D  fields  is  that  the  differences  between  the  passive  and  coupled  density 
fields  due  to  the  shift  in  the  equilibrium  point  is  eliminated,  since  this  shift  is  primarily 
horizontal.  These  profiles  also  provide  more  information  than  the  time  series  of  the  extremum 
density  values.  Such  profiles  for  passive,  coupled,  and  the  difference  fields  for  all  bubble  sizes 
are  shown  in  Fig.  6.  Here  we  observe  that  as  the  bubble  size  decre£ises  the  spike  representing 
accumulation  of  bubbles  decreases  ^nd  becomes  broader.  More  striking  is  reduction  in  the 
difference  profile  with  decreasing  bubble  size. 

6.2  !F  —  2:  Weak  gravity 

If  we  now  consider  larger  values  of  !F,  we  reduce  the  gravitational  effects  on  the  flow  and 
bubbles.  In  terms  of  the  bubble  parameters,  for  the  same  bubble  sizes  as  in  the  previous 
section  we  maintain  the  same  values  of  A,  but  have  smaller  values  of  W.  The  effect  these 
smaller  rise  velocities  have  on  the  accumulation  of  bubbles  in  the  passive  case  is  to  move 
the  equilibrium  points  closer  to  the  vortex  center  and  to  increase  the  entrapment  rate,  as 
seen  in  Fig.  7.  The  increase  in  entrapment  rate  is  most  noticeable  for  the  larger  bubbles, 
with  A  =  5.8,  where  the  maximum  number  density  reaches  values  of  ~  20^.  These  large 
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accumulation  rates  can  be  attributed  to  the  smaller  correction  to  the  Stokes  drag  due  to  the 
proximity  of  the  equilibrium  points  to  the  vortex  center. 

In  the  coupled  simulations  with  =  2,  we  also  observe  a  larger  accumulation  of  bubbles 
than  in  the  coupled  simulations  with  ^  =  1.  One  might  have  expected  that  the  large  accu¬ 
mulations,  especially  in  the  A  =  5.8  case,  would  alter  the  pressure  field  to  the  extent  that  the 
accumulation  rate  would  be  drastically  reduced,  even  to  the  point  where  the  accumulation 
itself  decreases,  as  in  the  latter  stages  in  the  =  1  and  A  =  5.8  case.  However,  when 
comp^ing  simulations  with  different  Froude  numbers  we  must  remember  that  the  effect  of 
bubble  accumulation  on  the  pressure  and  vorticity  fields  is  dependent  not  only  on  the  mag¬ 
nitude  of  the  density  gradients,  but  also  in  the  Froude  number  through  the  coefficient. 
As  a  result,  although  we  see  a  greater  accumulation  of  bubbles  for  the  T  =  1  case  than  in 
the  'J-  —  \  case,  the  pressure  graxiients  for  the  ^  =  2,«4  =  5.8  case,  shown  in  Fig.  8,  are 
aiffected  less  by  the  bubble  accumulation  than  in  the  =  1  case  of  Fig.  2b.  From  Fig.  8  we 
notice  that  the  change  in  the  vorticity  field  between  passive  and  coupled  cases  is  larger  in 
magnitude  than  for  the  J-  =\  case,  although  the  difference  is  confined  to  a  smaller  region. 
Therefore  we  have  a  lajge  difference  in  vorticity  near  the  vortex  center,  but  elsewhere  the 
vorticity  field  is  not  affected. 

The  compactness  of  the  accumulation  region  can  be  attributed  to  several  different  reasons. 
The  most  obvious  is  due  to  the  increase  in  the  entrapment  rate,  once  again  resulting  from 
the  small  slip  velocity  and  hence  drag  correction  when  the  equilibrium  point  lies  near  the 
vortex  center.  The  second  reason  concerns  the  nature  of  the  equilibrium  point.  Up  until 
now  we  have  used  the  term  “equilibrium  point”  rather  loosely.  In  the  present  case,  where 
the  accumulation  occurs  near  the  vortex  center,  the  flow  is  well  approximated  by  the  solid- 
body  vortex  model,  and  as  a  result  the  equilibrium  point  is  a  single  point.  As  the  bubble 
accumulation  moves  further  away  from  the  vortex  center,  the  flow  less  resembles  the  solid- 
body  vortex  model,  and  eis  a  result  we  can  not  guarantee  the  uniqueness  of  an  equilibrium 
point,  i.e.  there  can  be  an  equilibrium  region  instead  of  an  equilibrium  point. 

As  we  look  at  smaller  bubbles  for  the  T  —2  case,  the  the  accumulation  rate  decreaises 
for  both  passive  and  coupled  c<ises.  Although  tiwax  is  greater  here  than  in  their  T  =  \ 
counterparts,  the  smaller  coefficient  of  the  density  gradient  terms  in  the  vorticity  and  pressure 
equations  begins  to  take  effect,  resulting  in  very  small  differences  between  the  passive  and 
coupled  cases,  as  seen  in  Fig.  9.  In  addition  to  the  small  differences  in  the  vorticity  and 
pressure  fields,  we  also  observe  simil<ir  trends  in  the  density  fields  between  passive  and 
couples  cases.  In  fact,  the  large  difference  in  the  density  field  between  the  passive  and 
coupled  cases  in  Fig.  8  are  attributed  primarily  to  the  slight  shift  in  the  equilibrium  point. 
We  can  get  a  better  view  of  this  by  exzmaining  the  vertical  density  profiles  in  Fig.  10.  Here 
we  see  that  although  the  density  near  the  equilibrium  points  is  much  smaller  th<m  in  the 

—  \  case  (note  the  different  horizontal  scale  in  Fig.  6  and  Fig.  10),  the  difference  between 
the  pcissive  and  coupled  simulations  for  edl  bubbles  sizes  is  much  smaller  for  7-  =2. 

To  summarize,  when  increasing  the  Froude  number  we  observe  a  greater  accumulation 
of  bubbles  resulting  from  the  reduced  drag  correction  near  the  equilibrium  points  and  the 
smellier  rise  velocities  which  allow  bubbles  more  time  to  become  entrapped  by  the  vortex. 
Although  we  see  large  density  gradients  in  both  passive  and  coupled  cases,  the  effect  on  the 
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vorticity  and  pressure  is  not  as  large  as  in  the  ^  =  1  simulations  due  to  the  coefficient 
of  the  density  gradient  terms  in  the  vorticity  and  presstire  equations,  with  the  only  exception 
being  the  vorticity  field  of  the  A  =  5.8  case. 


6.3  T  —  0.707:  Strong  gravity 

One  might  expect  from  the  results  with  T  —  "1  and  1  that  as  we  consider  T  <\,  where  grav- 
itationad  effects  are  larger,  we  would  observe  smaller  accumulation  rates  but  with  enhanced 
vorticity  production.  In  the  previous  sections  the  coefficient  of  the  vorticity  production  term, 
wais  for  the  most  part  dominant  over  the  changes  in  the  accumulation  rate  of  bubbles. 
Although  we  do  see  the  accumulation  of  bubbles  decreasing  with  increasing  rise  velocity, 
we  observe  another  phenomenon  which  reduces  vorticity  production.  This  is  depicted  in 
Figs.  ll(a)-(c).  At  time  f  =  30  we  observe  the  accumulation  of  bubbles  similar  to  the  cases 
with  large  Froude  numbers,  but  at  times  t  =  40,50  the  bubbles  leave  this  accumulation 
region  as  they  advect  with  the  fluid  and  rise  velocities.  The  pressure  gradients  near  the 
vortex  center  are  smaller  at  these  later  times,  and  as  a  result  the  vortex  is  no  longer  able 
to  trap  bubbles.  Note  that  this  occurs  in  both  passive  and  coupled  cases,  and  therefore  is 
not  a  result  of  the  two-way  coupling  between  phases.  A  similar  phenomenon  was  observed 
in  passive  bubbles  rising  through  steady-state  Stuau^t  vortices.^®  In  the  Stuairt  vortex  flow 
for  cases  without  equilibrium  points  there  exists  a  point  which  is  least  unstable  amd  allows 
bubbles  to  be  retained  by  the  flow,  but  not  trapped.  The  motion  of  bubbles  in  our  unsteady 
flow  results  from  both  the  disappeau-ance  of  the  equilibrium  point  and  to  a  lesser  extent  the 
retention  of  bubbles  near  a  least  unstable  point.  As  the  bubbles  axe  released  from  this  region 
they  form  streaks  which  advect  with  the  fluid  and  rise  velocities. 

The  effect  of  these  streaks  on  the  maximum  bubble  concentration  is  shown  in  Fig.  12. 
For  the  largest  bubbles  the  maximixm  number  density  only  reaches  ~  1.5q.  After  t  30  the 
maximum  number  density  actuailly  decreaises  for  both  paissive  and  coupled  cases.  Because  of 
these  small  accumulations,  we  see  from  Figs.  ll(a)-(c)  that  there  is  little  difference  between 
the  vorticity  and  pressure  fields  of  the  passive  and  coupled  simulations.  As  a  result,  the 
number  density  fields  of  both  passive  and  coupled  simulations  are  roughly  the  same.  Note 
that  unlike  cases  with  stable  equilibrium  points,  the  vorticity  production  in  this  case  occurs 
along  the  streaks  of  bubbles  md  not  near  the  vortex  center. 

As  we  consider  smaller  bubbles  at  the  same  Froude  number,  the  decrease  in  rise  velocity 
accompcinying  the  difference  in  size  results  in  the  reappearance  of  equilibrium  point  during 
the  latter  stages.  As  a  result  we  observe  larger  accumulations  of  bubbles,  with  greater 
differences  between  passive  and  coupled  cases,  as  demonstrated  in  Figs.  12  and  13.  However, 
once  the  rise  velocity  is  small  enough  to  create  an  equilibrium  point,  as  in  the  case  for 
A  =  10.3,  then  as  we  move  to  smaller  bubbles  and  thus  rise  velocities,  such  as  those  with 
A  =  20,  we  once  ag2un  find  that  the  accumulation  decreases  due  to  lanrger  drag,  and  the 
difference  between  passive  and  coupled  cases  becomes  smaller.  This  suggests  that  at  some 
intermediate  bubble  size  there  is  an  optimum  coupling  effect.  For  large  bubbles,  there  is  no 
.equilibrium  point  and  hence  no  large  accumulation,  and  for  small  bubbles  the  accumulation 
rate  into  the  equilibrium  point  is  so  smadl  that  large  density  gradients  axe  not  produced.  The 
reeison  no  such  trend  was  observed  for  larger  Froude  numbers  is  due  to  restriction  on  bubble 
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size  needed  to  satisfy  our  approximations. 


7  Discussion  and  Conclusions 

Through  direct  numerical  simulations  we  have  obtained  information  regarding  the  effect 
of  two-way  coupling  of  bubbly  flows  with  dilute  concentrations  of  bubbles.  In  general  we 
are  interested  in  how  the  two-way  coupled  simulations  affect  both  bubble  and  flow  motion 
relative  to  the  passive  case.  To  summarize  the  results  of  the  previous  sections,  we  find  that 
the  vorticity  near  the  center  of  the  vortex  is  reduced  as  a  result  of  the  accumulation  of 
bubbles,  and  that  the  pressure  gradients  near  the  vortex  center  are  also  reduced,  resulting 
in  a  smaller  accumulation  for  the  coupled  simulations  relative  to  the  passive  case.  We 
have  explored  how  the  bubble  size  and  Froude  number  affect  these  changes.  In  general,  we 
notice  smaller  coupling  effects  for  smaller  bubbles,  since  the  accumulation  rate  decreases 
with  bubble  size.  For  most  of  the  the  simulations  considered  here,  the  effect  of  the  Froude 
number  on  the  results  is  dominated  by  the  coefficient  in  the  vorticity  production  and 
pressure  equations  rather  than  the  effect  of  the  Froude  number  on  the  accumulation  rate. 

Although  we  have  only  considered  the  bubble  size  and  Froude  numbers  as  parameters 
in  this  study,  the  information  we  obtained  from  this  can  explain  how  changes  in  other 
parameters  affect  the  results.  We  have  used  a  constant  void  fraction  in  this  study,  but  if  we 
substitute  the  equation  for  the  perturbation  density  into  the  momentum  equation,  we  see 
that  the  coefficient  of  the  coupling  terms  can  be  expressed  as  Thus  changes  in  the 

global  void  fraction  can  be  explained  directly  from  the  previous  results.  Moderate  changes 
in  the  Reynolds  number  to  not  effect  the  flow  greatly,  but  would  alter  the  values  of  A  and 
W  via  Eqs.  (12)  md  13. 

In  addition  to  the  effects  of  different  parameters  on  the  simulations,  we  should  also 
consider  how  different  configurations  may  affect  the  results.  In  this  study  we  have  only 
presented  material  corresponding  to  flows  uniformly  seeded  with  uniform  bubble  sizes.  We 
have  also  run  simulations  with  variable  bubbles  sizes,  ranging  the  entire  scope  of  bubbles 
considered  in  this  study.  In  general,  the  results  are  similar  to  what  is  observed  for  uniform 
bubbles  sizes,  except  that  the  accumulation  region  is  not  as  compact  and  the  changes  in  the 
pressure  and  vorticity  fields  not  as  great  due  to  the  smaller  density  gradients. 

Regardless  of  the  values  of  the  parameters,  one  recurring  theme  in  the  simulations  is 
that  in  most  cases  an  equilibrium  point  exists  in  the  flow  where  bubbles  accumulate.  The 
presence  of  sucij  a  point  might  seem  cc' trary  to  the  assumption  of  a  dilute  bubbly  flow.  We 
have  chosen  a  small  void  fraction  and  small  enough  bubbles  in  hopes  of  avoiding  violation  of 
the  dilute  assumption,  and  we  now  give  quantitative  justification  that  this  is  so.  A  violation 
of  the  weakly-dilute  assiimpcion  would  imply  that  bubbles  are  close  enough  that  bubbh  - 
bubble  interactions  become  important.  Although  we  know  of  no  precise  minimum  cutoff  for 
the  bubble  separation  to  exclude  any  such  interaction,  we  cam  easily  calculate  the  minimal 
bubble  separation  occurring  in  the  flow  at  any  time  simply  from  knowledge  of  the  maximum 
void  fraction  or  number  density.  The  separation  between  bubble  centers.  As,  is  given  by: 


At  the  begiiuiing  of  the  simulations,  with  e  s  1,  we  have  a  unifonn  separation  of  As/a  ~  7.4. 
For  all  but  one  simulation  the  maximum  void  fraction  is  ~  5e  giving  As/a  s  4.3.  The  only 
exception  is  the  case  with  T  =  2  and  A  =  5.8,  which  results  in  e  ~  20c  or  As/a  =  2.7.  For 
this  latter  case  the  dilute  assumption  is  most  certainly  violated  at  the  equilibrium  point,  but 
in  all  other  cases  the  minimum  separation  is  large  enough  for  the  bubble-bubble  interactions 
to  be  neglected  as  a  first  approximation. 
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Table  I.  Flow  and  bubble  parameters. 
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II 
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10.3 
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294,912 
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294,912 
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1  Time  sequence  of  the  vorticity  field  for  passive  cases.  Contour  levels  are  in 
increments  of  0.1,  where  solid  lines  represent  positive  and  dashed  lines  negative 
values  of  vorticity. 

2  Vorticity,  density  and  pressure  fields  for  the  coupled  and  passive  cases,  along 
with  their  differences,  at  a)  f  =  30  and  b)  t  =  50  for  =  1  and  A  =  5.8. 
The  solid  lines  indicate  positive  values  and  the  dashed  line  negative  values. 
The  increment  between  countour  levels  is  given  in  each  figure.  At  t  =  30, 
the  accumulation  is  not  large  enough  to  produce  substantial  changes  from  the 
coupling,  however  at  t  =  50  we  observe  differences  in  all  fields,  most  notably 
the  reduction  in  vorticity,  pressure  gradient,  and  bubble  accumulation  near 
the  vortex  center. 

3  PDF  of  relative  vorticity  between  passive  and  coupled  cases  at  t  =  50  for 
^  =  1  and  A  =  5.8.  The  small  peak  at  laxge  negative  values  indicates  the 
intensity  of  the  vorticity  difference  near  the  vortex  center. 

4  Vorticity,  density  and  pressure  fields  for  the  coupled  and  passive  cases,  along 
with  their  differences,  at  t  =  50  for  .^  =  1  and  A  =  10.3.  For  these  smaller 
bubbles  we  observe  reduced  accumulation  and  hence  weaker  coupling. 

5  Time  series  of  maximum  and  minimum  number  density,  77,  for:  A  =  5.8 
(top),  10.3  (center),  and  23  (bottom).  All  cases  have  ^  =s  1.  The  solid  lines 
represent  two-way  coupled  and  dashed  lines  passive  simulations.  These  times 
series  show  that  the  accumulation  of  bubbles  decreases  with  bubble  size(larger 
.A).  Note  that  for  A  =  5.8  the  maximum  accumulation  for  the  coupled  case 
becomes  smaller  as  bubbles  are  released  &om  the  equilibrium  points. 

6  Vertical  density  profiles  averaged  over  horizontal  slices  at  t  =  50  for:  A  — 
5.8  (top),  10.3  (center),  and  23  (bottom).  All  cases  have  ^  =  1.  Here  we 
observe  how  the  acciimulation  and  difference  between  coupled  and  passive 
cases  decrease  with  bubble  size. 

7  Time  series  of  maximum  and  minimum  number  density  for:  .<4  =  5.8  (top), 
10.3  (center),  and  23  (bottom).  All  cases  have  T  —  1.  The  solid  lines  rep¬ 
resent  two-way  coupled  and  dashed  lines  passive  simulations.  The  reduced 
gravitationcd  effects  result  in  much  larger  accumulations  for  both  passive  and 
coupled  cases  than  for  the  T  —  \  cases. 

8  Vorticity,  density  and  pressure  fields  for  the  coupled  and  passive  cases,  along 
with  their  differences,  at  t  =  50  for  .F  =  2  and  A  —  5.8.  For  this  larger  Froude 
number  emd  hence  smaller  rise  velocity,  the  bubbles  accumulate  more  towards 
the  vortex  center  with  an  enhanced  acciimulation  rate. 

9  Vorticity,  density  and  pressure  fields  for  the  coupled  and  passive  cases,  along 
with  their  differences,  at  t  =  50  for  =  2  and  A  =  10.3.  Although  the 
accumulation  is  greater  than  for  the  saune  sized  bubbles  with  7  —  \  due  to 
the  coeflBicient  the  coupling  effects  are  smaller. 
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10  Vertical  density  profiles  averaged  over  horizontal  slices  at  t  =  50  for:  A  =  5.8 
(top),  10.3  (center),  and  23  (bottom).  All  cases  have  T  =  2.  These  Passive- 
Coupled  column  clearly  shows  the  weaker  coupling  effects  relative  to  Fig.  6. 

11  Vortidty,  density  and  pressure  fields  for  the  coupled  and  passive  cases,  along 
with  their  differences,  for  A  —  5.8  and  T  =  0.707  at  (a)  f  =  30,  (b)  i  =  40,  (c) 
t  =  50.  Although  we  see  the  accumulation  of  bubbles  similar  to  previous  cases 
at  t  =  30,  at  later  times  the  pressure  gradient  is  not  large  enough  to  entrap 
bubbles,  and  as  a  result  we  obserrve  streaks  of  bubbles  which  are  produced 
from  a  least  unstable  point  replacing  the  equilibrium  point. 

12  Time  series  of  maximum  and  minimum  number  density  for:  A  =  5.8  (top), 
10.3  (center),  and  23  (bottom).  AU  cases  have  ^  =  0.707.  The  solid  lines 
represent  two-way  coupled  and  dashed  lines  passive  simulations.  The  larger 
gravitational  effects  result  in  smaller  accumulations,  most  noticeable  in  the 
j4,  =  5.8  case  where  there  is  no  equilibrium  point  for  the  bubbles  at  later  times. 

13  Vertical  density  profiles  averaged  over  horizontal  slices  at  t  =  50  for:  A  =  5.8 
(top),  10.3  (center),  and  23  (bottom).  All  cases  have  =  .707.  This  suggests 
that  there  is  an  optimum  bubble  size  in  terms  of  bubble  accumulation  for 
cases  with  small  For  the  smallest  bubbles  (bottom)  the  accumulation  rate 
is  small,  but  for  large  bubbles  (top)  rise  velocity  is  so  large  that  equilibrium 
points  do  not  exist. 
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Figure  1:  Time  sequence  of  the  vortidty  field  for  passive  cases.  Contotir  levels  are  in  incre¬ 
ments  of  0.1,  where  solid  lines  represent  positive  and  dashed  lines  negative  values  of  vorticity. 
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Figure  2:  Vorticity,  density  and  pressure  fields  for  the  coupled  and  passive  cases,  along  with 
their  differences,  at  a)  f  =  30  and  b)  t  =  50  for  ^  =  1  and  A  =  5.8.  The  solid  lines  indicate 
positive  values  and  the  dashed  line  negative  values.  The  increment  between  countour  levels  is 
given  in  each  figure,  kit  =  30,  the  accumulation  is  not  large  enough  to  produce  substantial 
changes  from  the  coupling,  however  at  t  =  50  we  observe  differences  in  all  fields,  most  notably 
the  reduction  in  vorticity,  pressure  gradient,  and  bubble  accxunxilation  near  the  vortex  center. 
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Figure  3:  PDF  of  relative  vorticity  between  passive  and  coupled  cases  at  t  =  50  for  ^  =  1 
and  A  =  5.8.  Tbe  small  peak  at  large  negative  values  indicates  the  intensity  of  the  vorticity 
difference  near  the  vortex  center. 


Figure  4:  Vorticity,  density  and  pressure  fields  for  the  coupled  and  passive  cases,  along  with 
their  differences,  at  t  =  50  for  ^  =  1  and  A  =  10.3.  For  these  smaller  bubbles  we  observe 
reduced  accumulation  and  hence  weaker  coupling. 


Figure  5:  Time  series  of  maximum  and  minimum  number  density,  r},  for:  A  =  5.8  (top),  10.3 
(center),  and  23  (bottom).  All  cases  have  ^  =1.  The  solid  lines  represent  two-way  coupled 
and  dashed  lines  passive  simulations.  These  times  series  show  that  the  accumulation  of  bub¬ 
bles  decreases  with  bubble  size(larger  A).  Note  that  for  A  =  5.8  the  maximum  accumulation 
for  the  coupled  case  becomes  smaller  as  bubbles  are  released  from  the  equilibrium  points. 


Figure  6:  Vertical  density  profiles  averaged  over  horizontal  slices  at  t  =  50  for:  A  =  5.8 
(top),  10.3  (center),  and  23  (bottom).  All  cases  have  ^  =  1.  Here  we  observe  how  the 
accumulation  and  difference  between  coupled  emd  passive  cases  decrease  with  bubble  size. 
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Figure  7:  Time  series  of  maximum  and  minimum  number  density  for:  A  =  5.8  (top),  10.3 
(center),  and  23  (bottom).  All  cases  have  ^  =  2.  The  solid  lines  represent  two-way  coupled 
and  dashed  lines  passive  simulations.  The  reduced  gravitational  effects  result  in  much  larger 
accumulations  for  both  passive  and  coupled  cases  than  for  the  ^  =  1  rases. 
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Figure  8:  Vorticity,  density  and  pressure  fields  for  the  coupled  and  passive  cases,  along  with 
their  differences,  at  t  =  50  for  ^  =  2  and  A  =  5.8.  For  this  larger  Froude  number  and 
hence  smaller  rise  velocity,  the  bubbles  accumiilate  more  towards  the  vortex  center  with  an 
enhanced  accumulation  rate. 
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Figtire  9:  Vorticity,  density  and  pressure  fields  for  the  coupled  and  passive  cases,  along  with 
their  differences,  at  t  =  50  for  ^  =  2  and  A  —  10.3.  Although  the  accumulation  is  greater 
than  for  the  same  sized  bubbles  with  ^  =  1,  due  to  the  coefficient  the  coupling  effects 
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Figure  10:  Vertical  density  profiles  averaged  over  horizontal  slices  at  t  =  50  for  A  =  5.8 
(top),  10.3  (center),  and  23  (bottom).  All  cases  have  ^  =  2.  These  Passive-Coupled  column 
clearly  shows  the  weaker  coupling  effects  rdative  to  Fig.  6. 
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Figure  11:  Vorticity,  density  and  pressure  fields  for  the  coupled  and  passive  cases,  along  with 
their  differences,  for  =  5.8  and  ^  =  0.707  at  (a)  i  =  30,  (b)  t  =  40,  (c)  t  =  50.  Although 
we  see  the  accumulation  of  bubbles  similar  to  previous  cases  at  t  =  30,  at  later  times  the 
pressure  gradient  is  not  large  enough  to  entrap  bubbles,  and  as  a  result  we  obserrve  streaks 
of  bubbles  which  are  produced  from  a  least  unstable  point  replacing  the  equilibrium  point. 
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Figure  12:  Time  series  of  maximum  and  minimum  number  density  for  A  =  5.8  (top), 
10.3  (center),  and  23  (bottom).  All  cases  have  ^  =  0.707.  The  solid  lines  represent  two-way 
coupled  and  dashed  lines  pMsive  simulations.  The  larger  gravitational  effects  result  in  smaller 
accumulations,  most  noticeable  in  the  A  =  5.8  case  where  there  is  no  equilibrium  point  for 
the  bubbles  at  later  times. 
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Figure  13:  Vertical  density  profiles  averaged  over  horizontal  slices  at  t  =  50  for:  A  =  5.8 
(top),  10.3  (center),  and  23  (bottom).  All  cases  have  ^  =  .707.  This  suggests  that  there  is 
an  optimum  bubble  size  in  terms  of  bubble  accumulation  for  cases  with  small  For  the 
smallest  bubbles  (bottom)  the  accumulation  rate  is  small,  but  for  large  bubbles  (top)  rise 
velocity  is  so  large  that  equilibrium  points  do  not  exist. 
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